function [quality, quality_map] = img_qi(img1, img2, imgf, block_size)

%========================================================================
%
%Copyright (c) 2001 The University of Texas at Austin
%All Rights Reserved.
% 
%This program is free software; you can redistribute it and/or modify
%it under the terms of the GNU General Public License as published by
%the Free Software Foundation; either version 2 of the License, or
%(at your option) any later version.
% 
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
% 
%The GNU Public License is available in the file LICENSE, or you
%can write to the Free Software Foundation, Inc., 59 Temple Place -
%Suite 330, Boston, MA 02111-1307, USA, or you can find it on the
%World Wide Web at http://www.fsf.org.
%
%Author  : Zhou Wang 
%Version : 1.0
% 
%The authors are with the Laboratory for Image and Video Engineering
%(LIVE), Department of Electrical and Computer Engineering, The
%University of Texas at Austin, Austin, TX.
%
%Kindly report any suggestions or corrections to zwang@ece.utexas.edu
%
%Acknowledgement:
%The author would like to thank Mr. Umesh Rajashekar, the Matlab master
%in our lab, for spending his precious time and giving his kind help
%on writing this program. Without his help, this program would not
%achieve its current efficiency.
%
%========================================================================
%
%This is an efficient implementation of the algorithm for calculating
%the universal image quality index proposed by Zhou Wang and Alan C. 
%Bovik. Please refer to the paper: Zhou Wang and Alan C. Bovik "A Universal
%Image Quality Index" IEEE Signal Processing Letters, 2001. In order to
%run this function, you must have Matlab's Image Processing Toobox.
%
%You can change this program as you like and use it anywhere, but please
%refer to its original source (cite our paper and our web page at
%http://anchovy.ece.utexas.edu/~zwang/research/quality_index/demo.html).
%
%Input : an original image and a test image of the same size
%Output: (1) an overall quality index of the test image, with a value
%            range of [-1, 1].
%        (2) a quality map of the test image. The map has a smaller
%            size than the input images. The actual size is
%            img_size - BLOCK_SIZE + 1.
%
%Usage:
%
%1. Load the original and the test images into two matrices
%   (say img1 and img2)
%
%2. Run this function in one of the two ways:
%
%   % Choice 1 (suggested):
%   [qi qi_map] = img_qi(img1, img2);
%
%   % Choice 2:
%   [qi qi_map] = img_qi(img1, img2, BLOCK_SIZE);
%
%   The default BLOCK_SIZE is 8 (Choice 1). Otherwise, you can specify
%   it by yourself (Choice 2).
%
%3. See the results:
%
%   qi                    %Gives the over quality index.
%   imshow((qi_map+1)/2)  %Shows the quality map as an image.
%
%========================================================================

if (nargin == 1 | nargin > 4)
   quality = -Inf;
   quality_map = -1*ones(size(img1));
   return;
end

if ((size(img1) ~= size(img2)) | (size(img1) ~= size (imgf)))
   quality = -Inf;
   quality_map = -1*ones(size(img1));
   return;
end

if (nargin == 3)
   block_size = 8;
end

N = block_size.^2;
sum2_filter = ones(block_size);

img1_sq   = img1.*img1;
img2_sq   = img2.*img2;
imgf_sq   = imgf.*imgf;

img1f = img1.*imgf;
img2f = img2.*imgf;

img1_sum   = filter2(sum2_filter, img1, 'valid');
img2_sum   = filter2(sum2_filter, img2, 'valid');
imgf_sum   = filter2(sum2_filter, imgf, 'valid');

img1_sq_sum = filter2(sum2_filter, img1_sq, 'valid');
img2_sq_sum = filter2(sum2_filter, img2_sq, 'valid');
imgf_sq_sum = filter2(sum2_filter, imgf_sq, 'valid');

img1f_sum = filter2(sum2_filter, img1f, 'valid');
img2f_sum = filter2(sum2_filter, img2f, 'valid');

img1f_sum_mul = img1_sum.*imgf_sum;
img2f_sum_mul = img2_sum.*imgf_sum;

img1f_sq_sum_mul = img1_sum.*img1_sum + imgf_sum.*imgf_sum;
img2f_sq_sum_mul = imgf_sum.*imgf_sum + img2_sum.*img2_sum;

numerator1 = 4*(N*img1f_sum - img1f_sum_mul).*img1f_sum_mul;
denominator1 = N*(img1_sq_sum + imgf_sq_sum) - img1f_sq_sum_mul;
denominator1a = denominator1.*img1f_sq_sum_mul;

numerator2 = 4*(N*img2f_sum - img2f_sum_mul).*img2f_sum_mul;
denominator2 = N*(imgf_sq_sum + img2_sq_sum) - img2f_sq_sum_mul;
denominator2a = denominator2.*img2f_sq_sum_mul;

sf1=spfreqmat(img1,block_size);
sf2=spfreqmat(img2,block_size);
lamda=sf1./(sf1+sf2);
lamda(isnan(lamda))=0.5;

quality_map1 = ones(size(denominator1));
quality_map2 = ones(size(denominator1));

index1 = (denominator1 == 0) & (img1f_sq_sum_mul ~= 0);
index2 = (denominator2 == 0) & (img2f_sq_sum_mul ~= 0);

quality_map1(index1) = 2*img1f_sum_mul(index1)./img1f_sq_sum_mul(index1);
quality_map2(index2) = 2*img2f_sum_mul(index2)./img2f_sq_sum_mul(index2);

index1 = (denominator1a ~= 0);
index2 = (denominator2a ~= 0);

quality_map1(index1) = numerator1(index1)./denominator1a(index1);
quality_map2(index2) = numerator2(index2)./denominator2a(index2);

quality_map = lamda.*quality_map1 + (1-lamda).*quality_map2;

quality = mean2(quality_map);
